library(Seurat)
library(ggpubr)
library(viridis)
library(tidyr)
library(dplyr)
library(tidyr)
library(readxl)
library(stringr)
library(gridExtra)
library(grid)
library(purrr)
library(tidyverse)
out_dir = "/fh/fast/sun_w/zhaoheng_li/MSK_DEC24/Analysis/outs/differential_expression"
celllevel_out_dir = file.path(out_dir,"cell-level")
pseudobulk_out_dir = out_dir
fig_dir =  file.path(out_dir,"compare-DE")
for(tp in c("D-1","D3","D7","D11","D18")){
  print(paste0("time point: ", tp))
  pseudobulk_res = read.delim(file.path(pseudobulk_out_dir, paste0("DE_pseudobulk_", tp, ".tsv")), 
                   header = TRUE, sep = "\t", stringsAsFactors = FALSE)
  
  wilcox_res = read.csv(file.path(celllevel_out_dir,paste0("DE_cell-level_",tp,".csv")),row.names = 1)
  merged <- pseudobulk_res %>%
    select(celltype, genotype, gene, pvalue) %>%
    inner_join(
      wilcox_res %>% select(celltype, genotype, gene, p_val),
      by = c("celltype", "genotype", "gene")
    ) %>%
    mutate(
      pvalue   = pmax(pvalue, .Machine$double.xmin),
      p_val    = pmax(p_val,  .Machine$double.xmin),
      nl10_pb  = -log10(pvalue),
      nl10_wc  = -log10(p_val)
    )
  for(ct in unique(merged$celltype)){
    print(paste0("cell type: ", ct))
    fig = merged %>%
      filter(celltype == ct) %>%
      ggplot(aes(x = nl10_pb, y = nl10_wc)) +
      geom_point(alpha = 0.6, size = 1) +
      geom_abline(slope = 1, intercept = 0, linetype = "dashed") +
      stat_cor(method = "spearman", label.x = 1, label.y = 19, size = 3) + 
      coord_cartesian(xlim = c(0, 20), ylim = c(0, 20)) +
      labs(
        x = "-log10(pseudobulk p-val)",
        y = "-log10(wilcoxon p-val)",
        title = paste(tp, ct, sep = "_")
      ) +
      theme_classic() +
      facet_wrap(~genotype)
    print(fig)
    ggsave(file.path(fig_dir,paste0(tp,"_",ct,".png",sep='')),fig,width = 12,height = 12)
  }
}
## [1] "time point: D-1"
## [1] "cell type: ESC"

## [1] "time point: D3"
## [1] "cell type: DE"

## [1] "cell type: ESC_DE"

## [1] "time point: D7"
## [1] "cell type: Liver"

## [1] "cell type: PFG"

## [1] "time point: D11"
## [1] "cell type: PP"

## [1] "time point: D18"
## [1] "cell type: Ductal"

## [1] "cell type: EnP"

## [1] "cell type: SC-alpha"

## [1] "cell type: SC-beta"

## [1] "cell type: SC-delta"

## [1] "cell type: SC-EC"

## [1] "cell type: Stromal"

Session information

gc()
##            used  (Mb) gc trigger  (Mb)  max used  (Mb)
## Ncells  5461632 291.7   11419752 609.9  11419752 609.9
## Vcells 42015379 320.6  108915709 831.0 108915709 831.0
sessionInfo()
## R version 4.4.0 (2024-04-24)
## Platform: x86_64-pc-linux-gnu
## Running under: Ubuntu 18.04.6 LTS
## 
## Matrix products: default
## BLAS/LAPACK: FlexiBLAS OPENBLAS;  LAPACK version 3.11.0
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## time zone: America/Los_Angeles
## tzcode source: system (glibc)
## 
## attached base packages:
## [1] grid      stats     graphics  grDevices utils     datasets  methods  
## [8] base     
## 
## other attached packages:
##  [1] lubridate_1.9.3    forcats_1.0.0      readr_2.1.5        tibble_3.2.1      
##  [5] tidyverse_2.0.0    purrr_1.0.2        gridExtra_2.3      stringr_1.5.1     
##  [9] readxl_1.4.3       dplyr_1.1.4        tidyr_1.3.1        viridis_0.6.5     
## [13] viridisLite_0.4.2  ggpubr_0.6.0       ggplot2_3.5.1      Seurat_5.1.0      
## [17] SeuratObject_5.0.2 sp_2.1-4          
## 
## loaded via a namespace (and not attached):
##   [1] RColorBrewer_1.1-3     rstudioapi_0.16.0      jsonlite_1.8.8        
##   [4] magrittr_2.0.3         spatstat.utils_3.0-4   farver_2.1.2          
##   [7] rmarkdown_2.26         ragg_1.3.1             vctrs_0.6.5           
##  [10] ROCR_1.0-11            spatstat.explore_3.2-7 rstatix_0.7.2         
##  [13] htmltools_0.5.8.1      broom_1.0.5            cellranger_1.1.0      
##  [16] sass_0.4.9             sctransform_0.4.1      parallelly_1.37.1     
##  [19] KernSmooth_2.23-22     bslib_0.7.0            htmlwidgets_1.6.4     
##  [22] ica_1.0-3              plyr_1.8.9             plotly_4.10.4         
##  [25] zoo_1.8-12             cachem_1.0.8           igraph_2.0.3          
##  [28] mime_0.12              lifecycle_1.0.4        pkgconfig_2.0.3       
##  [31] Matrix_1.7-0           R6_2.5.1               fastmap_1.2.0         
##  [34] fitdistrplus_1.2-1     future_1.33.2          shiny_1.8.1.1         
##  [37] digest_0.6.35          colorspace_2.1-0       patchwork_1.2.0       
##  [40] tensor_1.5             RSpectra_0.16-1        irlba_2.3.5.1         
##  [43] textshaping_0.3.7      labeling_0.4.3         progressr_0.14.0      
##  [46] timechange_0.3.0       fansi_1.0.6            spatstat.sparse_3.0-3 
##  [49] httr_1.4.7             polyclip_1.10-6        abind_1.4-5           
##  [52] compiler_4.4.0         withr_3.0.0            backports_1.4.1       
##  [55] carData_3.0-5          fastDummies_1.7.3      highr_0.10            
##  [58] ggsignif_0.6.4         MASS_7.3-60.2          tools_4.4.0           
##  [61] lmtest_0.9-40          httpuv_1.6.15          future.apply_1.11.2   
##  [64] goftest_1.2-3          glue_1.7.0             nlme_3.1-164          
##  [67] promises_1.3.0         Rtsne_0.17             cluster_2.1.6         
##  [70] reshape2_1.4.4         generics_0.1.3         gtable_0.3.5          
##  [73] spatstat.data_3.0-4    tzdb_0.4.0             hms_1.1.3             
##  [76] data.table_1.15.4      car_3.1-2              utf8_1.2.4            
##  [79] spatstat.geom_3.2-9    RcppAnnoy_0.0.22       ggrepel_0.9.5         
##  [82] RANN_2.6.1             pillar_1.9.0           spam_2.10-0           
##  [85] RcppHNSW_0.6.0         later_1.3.2            splines_4.4.0         
##  [88] lattice_0.22-6         survival_3.6-4         deldir_2.0-4          
##  [91] tidyselect_1.2.1       miniUI_0.1.1.1         pbapply_1.7-2         
##  [94] knitr_1.46             scattermore_1.2        xfun_0.44             
##  [97] matrixStats_1.3.0      stringi_1.8.4          lazyeval_0.2.2        
## [100] yaml_2.3.8             evaluate_0.23          codetools_0.2-20      
## [103] cli_3.6.2              uwot_0.2.2             systemfonts_1.1.0     
## [106] xtable_1.8-4           reticulate_1.36.1      munsell_0.5.1         
## [109] jquerylib_0.1.4        Rcpp_1.0.12            globals_0.16.3        
## [112] spatstat.random_3.2-3  png_0.1-8              parallel_4.4.0        
## [115] dotCall64_1.1-1        listenv_0.9.1          scales_1.3.0          
## [118] ggridges_0.5.6         leiden_0.4.3.1         rlang_1.1.3           
## [121] cowplot_1.1.3